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Abstract 

We analyze the existence and stability of localized solutions in the one-dimensional 
discrete nonlinear Schrodinger (DNLS) equation with a combination of competing 
self-focusing cubic and defocusing quintic onsite nonlinearities. We produce a stabil- 
ity diagram for different families of soliton solutions, that suggests the (co) existence 
of infinitely many branches of stable localized solutions. Bifurcations which occur 
with the increase of the coupling constant are studied in a numerical form. A varia- 
tional approximation is developed for accurate prediction of the most fundamental 
and next-order solitons together with their bifurcations. Salient properties of the 
model, which distinguish it from the well-known cubic DNLS equation, are the ex- 
istence of two different types of symmetric solitons and stable asymmetric soliton 
solutions that are found in narrow regions of the parameter space. The asymmet- 
ric solutions appear from and disappear back into the symmetric ones via loops of 
forward and backward pitchfork bifurcations. 
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1 Introduction 



Discrete nonlinear Schrodinger (DNLS) equations constitute an important 
class of discrete lattice models that are of great interest in their own right [1], 
and also find direct applications to the description of arrays of waveguides in 
nonlinear optics, as predicted in Ref . [2] , and for the first time realized exper- 
imentally in Ref. [3], which used a set of parallel semiconductor waveguides 
made on a common substrate (see review article Ref. [4] and further refer- 
ences therein). In optics, quasi-discrete waveguide arrays can also be created 
as virtual photonic lattices in photorefractive materials (see review Ref. [5] 
and references therein), and can be approximately described by the DNLS 
equations. The waveguide arrays may support both spatial solitons [3,4] and 
quasi-discrete spatiotemporal collapse [6]. 

Besides nonlinear optics, the DNLS model also describes a Bose-Einstein con- 
densate trapped in a strong optical lattice (sinusoidal potential acting on 
atoms in the condensate), as predicted theoretically [7] and observed in the 
experiment [8] (see also the recent review Ref. [9]). Additionally, DNLS equa- 
tions may be naturally derived, in the rotating-phase approximation, from 
various nonlinear-lattice models that give rise to discrete breathers (alias in- 
trinsic localized modes), see theoretical papers Refs. [10] and [11], and first 
reports of the experimental making of these breathers in Ref. [12]. 

Properties of discrete solitons in the DNLS equation with the simplest, cubic, 
nonlinearity have been studied in detail, including three dimensional settings 
[13], and are now well understood. These solitons were experimentally observed 
in arrays of nonlinear optical waveguides [3,4]. They also correspond, in the 
DNLS approximation, to the intrinsic localized modes in more sophisticated 
dynamical lattices. 

Nonlinear Schrodinger (NLS) equations with more complex nonlinearities were 
studied in detail in continuum models. As well as their cubic counterparts, such 
models are of interest by themselves, and may also have direct applications 
[14]. In particular, glasses and organic optical media whose dielectric response 
features the cubic-quintic (CQ) nonlinearity, i.e., a self-defocusing quintic cor- 
rection to the self-focusing cubic Kerr effect, are known [15]. Properties of 
solitons in the NLS equations with the CQ nonlinearity may be very different 
from those in the simplest cubic equation, especially in the case when the 
higher-order nonlinearity is combined with a periodic potential. Recently, a 
great variety of multistable solitons with different numbers of peaks and dif- 
ferent symmetries (even, odd, etc.) have been found in the CQ NLS equation 
embedded in the linear potential of the Kronig-Penney (KP) type (a periodic 
array of rectangular potential wells) [16], after bistable solitons were studied 
in the CQ NLS equation with a single rectangular potential well [17]. Solitons 
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in the continuum cubic NLS equation with the KP potential have been studied 
too [18]. 



The limiting case of the CQ NLS equation with a very strong KP potential 
naturally reduces to the DNLS equation with the CQ nonlinearity, and our ob- 
jective in this paper is to construct solitons in this discrete model and explore 
their stability. The model is not only interesting by itself (as is shown in the 
present work), but may also be realized experimentally in the form of an array 
of waveguides built of the above-mentioned optical materials featuring the CQ 
nonlinearity [15]. It is relevant to mention that stable discrete solitons were 
recently found in the DNLS equation with saturable nonlinearity [20,21] (note 
that the latter model was introduced back in 1975 by Vinetskii and Kukhtarev 
[22]), and, moreover, optical discrete solitons supported by the saturable self- 
defocusing nonlinearity were experimentally created using the photovoltaic 
effect in a waveguiding lattice built into a photorefractive crystal [23]. We 
show in this work that discrete solitons in the CQ model are very different 
from their counterparts investigated in the aforementioned works [20,21,23], 
(most importantly, they feature a great multistability, as shown below) due to 
the fundamental fact that, unlike the saturable nonlinearity, the combination 
of the CQ terms features competition of self-focusing and defocusing types. 
For the same reason, the solitons in the CQ DNLS equation are drastically 
different from ones investigated earlier [24,25] in the DNLS equation with a 
single nonlinear term of an arbitrary power (for instance, quintic instead of 
cubic). Finally, a quantum version of a finite-length DNLS equation (Bose- 
Hubbard model) with the CQ nonlinearity and periodic boundary conditions 
was considered in Ref. [26], where states with a small number q of quanta (in 
most cases, q < 6) were constructed. 

The manuscript is organized as follows. In the next section, we introduce the 
model and its stationary solutions, and derive a two-dimensional map to gen- 
erate discrete solitons corresponding to homoclinic solutions. Section 3 reports 
various (multistable) soliton solutions and their stability. In Section 4, we fo- 
cus on bifurcations that create/annihilate different solutions and account for 
exchange of stability between them. In Section 5 we present an analytical vari- 
ational approximation that correctly predicts the main bifurcation branches. 
Section 6 concludes the paper. 



2 The model and dynamical reductions 

The one-dimensional cubic-quintic discrete nonlinear Schrodinger (CQ DNLS) 
equation is 
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where ip n are the complex fields at site n (in the case of the above-mentioned 
array of optical waveguides, ip n is amplitude of the electromagnetic wave in 
the given core), ip = dip/alt (in the above-mentioned model of the waveguide 
array, the evolutional variable is actually not time but the coordinate along the 
waveguide), and the discrete second derivative (discrete-diffraction operator 
in the array of waveguides) is CA^>„ = C (ip n +i + V'n-i — 2ip n ), where C is the 
constant of the tunnel coupling between the cores. The third and fourth terms 
in Eq. (1) represent, respectively, the cubic and quintic nonlinearities. We 
assume C,B,Q > 0, which (as said above) corresponds to the most natural 
case of the self-focusing cubic (Kerr) nonlinearity competing with its self- 
defocusing quintic counterpart. Upon renormalization of ip and t, we set B = 2 
and Q = 1. 



Equation (1) conserves two dynamical invariants: the Hamiltonian, 



n 

and norm, 



CV>* (V>n+1 + i>n-l ~ 2^n) + l^nf ~ ^ \^nf 



(2) 



M = £ivg 2 (3) 



(in the application to the optical waveguide array, the latter is the total power 
of the light signal). 

We aim to look for a soliton solution with a frequency fi by substituting 

= u n exp(-ifit) (4) 

in Eq. (1). The real stationary lattice field u n must solve the equation 

fm n + C(u n+1 + u n -i - 2u n ) + 2u 3 n - u b n = 0, (5) 

supplemented by the condition of vanishing of u n at n — > ±oo, i.e., the soliton 
can be looked for as a homoclinic solution of Eq. (5) (an alternative approach 
would be to use an algebraic method of Ref. [27]). Note that stationary soliton 
solutions of Eq. (5) depend on two parameters, /i and C . 

We will study soliton solutions and their stability by viewing Eq. (5) as a 
recurrence relation between consecutive amplitudes, that can be cast in the 
form of a two-dimensional map, 

u n+l = au n -v n - 2C^ul + C~ l u\ 

(6) 

where a = 2 — fi/C. Constant solutions to Eq. (1) correspond to fixed points 
(FPs) of this map. There exists at most five FPs, that we arrange in increasing 
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Fig. 1. Stability diagram for the fixed points (FPs) of the two-dimensional map (6). 
In non shaded areas all the FPs are saddles. Shaded areas represent regions where 
the indicated FPs are centers. The diagonal line and the tilted parabola correspond, 
respectively, to C = /x/4 and (C + 1 + fi) 2 = 1 + \i. 

order and label as u^ + ,u , u , u ++ , with Uo = and 

u±± = ±\Jl±J\TJi. (7) 

The stability of all the FPs within the framework of map (6) can be easily 
derived from the linearization of the map around these FPs, leading to the 
stability chart displayed in Fig. 1. 

In this framework, soliton solutions correspond to homoclinic orbits connecting 
the FP at origin (u = 0) with itself, in the case when it is a saddle. The FP 
uq = is a saddle for {// < and C > /i/4} and for {/i > and C < /x/4}. 
We are only interested in physically meaningful couplings so we restrict our 
attention to C > 0. Furthermore, the regions for /i < — 1 and /i > produce 
stable and unstable manifolds that do not intersect each other. Therefore, in 
the remaining of this work, we restrict our area of interest to C > with 
-1< fi < 0. 

In Fig. 2, we depict a progression of the homoclinic tangles (webs of such or- 
bits) emanating from the saddle points as the coupling constant C increases. 
For small C, the homoclinic-connection structure is very rich, including many 
orbits (i.e., many soliton solutions). As C increases, many connections disap- 
pear through a series of bifurcations (see below), so that a single homoclinic 
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Fig. 2. Homoclinic tangles generating localized solutions of Eq. (1) for fi = —0.6, 
with the coupling constant C increasing from left to right and from top to bottom: 
(a) C = 0.15, (b) C = 0.4, (c) C = 0.8, and (d) C = 2. Saddle points and centers 
are designated by asterisks and circles, respectively. For small C, the homoclinic 
intersections have a much richer structure, including homoclinic and heteroclinic 
connections between all five fixed points. For large C, only one homoclinic and one 
heteroclinic solutions survive (together with their u n <-> —u n symmetric counter- 
parts). 

loop survives at very large C, which corresponds to the well-known exact 
soliton solution of the continuum CQ NLS equation [28]. 

Note that, in contrast with the cubic DNLS equation, the CQ discrete equa- 
tion gives rise to a pair of extra fixed points that support heteroclinic orbits 
(as shown in Fig. 2). The detailed study of these heteroclinic orbits, which 
correspond to dark solitons or kinks, falls outside the scope of this work. Here 
we concentrate on the homoclinic orbits and (bright) soliton solutions corre- 
sponding to them. 

In Fig. 3 we depict a selection of typical solitons corresponding to the homo- 
clinic tangles in Fig. 2. Each solution is numerically generated by taking the 
soliton shape as predicted in an approximate form by the homoclinic inter- 
sections, and then applying a Newton-type algorithm to find a numerically 
exact localized solution of Eq. (5). In panel (a) we show a couple of solitons 
coexisting at given parameter values. Panel (b) depicts a triplet of coexisting 
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Fig. 3. Typical examples of soliton solutions generated by the homoclinic tangles of 
Fig. 2 as the coupling parameter C increases. Insets pertain to the same values of 
parameters as the corresponding insets in Fig. 2. (a) Pair of coexisting symmetric 
solution for small C. (b) As C increases, only one stable symmetric solution exists 
(squares), together with a pair of asymmetric solutions (triangles and circles) ap- 
pearing in small regions of the parameter space (the asymmetric solutions have been 
vertically displaced for clarity of presentation), (c) Further increase in C destroys 
all stable solutions except for the symmetric one, that, as seen in panel (d), tends 
towards the continuum soliton in the limit of C — > oo. The parameters are fi = —0.6 
and (a) C = 0.15, (b) C = 0.4, (c) C = 0.8, (d) C = 2. 

solutions that form a part of a loop of pitchfork bifurcations responsible for 
the creation of asymmetric solutions (see below). Finally, panels (c) and (d) 
show the unique site-centered (the maximum of the soliton is located at a sin- 
gle central site, cf. Fig. 5(a)) solution which survives in the continuum limit, 
C -> oo. 



3 Multistability of discrete solitons 

We now aim to explore the structure of the homoclinic tangles and their bifur- 
cations in detail, varying the parameters /i and C . As previously mentioned, 
for small C the rich homoclinic structure leads to the coexistence of multiple 
solitons at the same values of fi and C, which is a distinctive feature of the 
CQ model with the competing nonlinearities: in the cubic DNLS equations, 
this variety of solitons is not observed [1]. Principal types of the localized so- 
lutions for small C (taking C = 0.1 as an example) are depicted in Fig. 4. 
The following scheme is adopted to denote different species of the solitons. 
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The solution generated by the first (main) homoclinic crossing of the stable 
and unstable manifolds of the origin (see Fig. 4(a)) is denoted by Si. This 
family corresponds to site-centered solitons, the label S standing for short, as 
the solitons of this type correspond to the shortest family of homoclinic cross- 
ings. As C is increased, the Si soliton suffers a series of alternating stability 
switches (bifurcations). We use the notation Sk to denote solitons correspond- 
ing to the series of stable regions between the stability switches. In this work 
we do not consider higher-order crossings corresponding to repeated iterates 
in both the stable and unstable manifolds, that would produce bound states of 
the discrete solitons, alias multi-humped or multibreather (non fundamental) 
solutions [29]. 




Fig. 4. Multistability of discrete solitons in the CQ DNLS equation is illustrated 
by displaying the homoclinic tangles generating stable localized solutions (see Fig. 
5). Parameters are [x = —0.6 and C = 0.1 (which corresponds to the point marked 
by the asterisk in the T\ region in Fig. 6). Panels (a), (b), (c) and (d) correspond, 
respectively, to the discrete solitons of the Si, T\, T%, and T3 types (see definitions 
in the text). 

The solitons generated by the homoclinic crossings in Figs. 4(b)- (d) corre- 
spond, in our notation, to T 1; T 2 and T 3 solitons, T standing for tall solitons. 
They are generated by the second family of crossings, and are characterized 
by a higher soliton maximum, in comparison with their S counterparts. 

The subscript in the notation for the Tk and Sk solutions also helps to dif- 
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ferentiate between site- centered solutions, for k odd, and bond-centered (alias 
inter site- centered) ones for even k, which feature two central sites with equal 
magnitude, see examples in Fig. 5. 




Fig. 5. Evolution of the solitons generated by the homoclinic intersections depicted 
in Fig. 4 after adding random perturbations, with a relative amplitude of 5%, to 
the initial configuration (dark line). 

The Si, Ti, T 2 and T 3 solitons generated by the highlighted crossings in Fig. 4 
are all stable solutions. The stability was checked by calculating the respective 
eigenvalues from Eq. (1) linearized around the stationary solutions, and also 
verified by direct numerical integration of the full equation (1) after adding a 
random perturbation to the soliton, with a (rather large) relative amplitude of 
5%. The evolution of the so perturbed solutions is displayed in Fig. 5, where 
the unperturbed solitons are shown by dark lines for t — 0. The perturbed 
solutions oscillates about the unperturbed solitons, confirming their stability 
(the oscillations do not fade because the system is conservative). 

In order to identify a stability region of the soliton family in the (/i, C) pa- 
rameter plane, we start with a particular soliton solution of the S or T type, 
generated as described above, and then continued it by varying the param- 
eters, simultaneously computing the stability eigenvalues. The continuation 
procedure started at small C where, as said above, it is easy to find a partic- 
ular solution. 

The resulting stability diagram for the solitons of the £1,2,3 and Ti j2) 3 types is 
displayed in Fig. 6. The Testable regions (shaded in Fig. 6) feature a tent-like 
shape, with the base on the line C = with — 1 < // < 0. Since all these regions 
share the common base, there should be a nontrivial area where, presumably, 
the Tk solitons exist and are stable for all values of k. It is also apparent from 
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Fig. 6. The stability diagram for single- hump localized solitons of the cubic-quintic 
discrete NLS equation (1). Note that each region has a "wedge" that penetrates 
into the corresponding region. The existence region of the solitons is bounded 
by the curve EB, at which the homoclinic connections of the origin disappear. The 
curves VAl and VA2 show the variational approximations described in Section 5. 
A more detailed description is given in the text. 

the stability diagram in Fig. 6 that the stability regions for the S& solutions 
with different k, in contrast to their counterparts, do not intersect each 
other. Therefore, the solitons of the S type feature no multistability. 

We note that each stability region features a wedge that penetrates into the 
region of stability of the Tj, solitons. Particularly, the 7\ region is completely 
embedded in the Si region. This property tends to suggest that, for any k, 
there always exists non-empty regions where the and S& solitons coexist 
and are simultaneously stable. It is interesting too that the S-stability regions 
outlive their T-counterparts as C increases. This is due to the fact that the T 
solutions correspond to the second family of homoclinic crossings that disap- 
pear, in saddle- node bifurcations (see below), earlier than the first family of 
the crossings (i.e., the S solutions). 

Existence regions for the solitons, which may be broader than the regions 
where the solitons are stable, are defined as those in which the stable and 
unstable manifolds emanating from the origin do intersect. We have computed 
such a region numerically. In Fig. 6, it is located to the right of the dashed- 
dotted line EB (existence boundary). To the left of this curve, no localized 
solutions are possible. As C increases, the EB curve approaches the line /x = 
—3/4, that corresponds to the existence border for solitons in the continuum 



10 





Rett.) 



0.5 




0.15 



Fig. 7. Evolution of unstable solitons (left panels), together with their instability 
spectrum (right panels). The top row depicts the evolution of an unstable S2 soliton, 
shown by the dashed line in Fig. 3(a) for (u,C) = (—0.6,0.15), initiated by a small 
perturbation. After some transient, the initial configuration (solid line) is changed 
by one featuring oscillations between two asymmetric states (dashed lines). The 
bottom row depicts the evolution of an unstable symmetric Si solution from Fig. 
3(b) for (//, C) = (—0.6,0.4). The evolution amounts to oscillations between the 
initial configuration (solid line) and its translation by one site (dashed line). The 
respective (in)stability spectra for the two solutions in the left column are depicted 
in the right column. 

limit of the CQ equation [16]. More specifically, the existence regions of the 
T solutions exactly coincide with their stability regions, i.e., the solitons of 
the T type are always stable. The existence region for the S solutions has 
a complex structure, contained in the region where homoclinic connections 
exists to the right of the EB curve, while the stability area for these solitons 
is really smaller than the existence region. We also note that (as explained in 
detail below) the Sk solutions represent only two distinct families of solitons, 
one site-centered (for odd k) and one bond-centered (for even k). 

Since some solitons of the S type are unstable, it is necessary to determine 
what the instability transforms them into. As an example, in the top row 
of Fig. 7 we display the evolution of an unstable S2 soliton in the stability 
region of T\ (for (//, C) = (—0.6,0.15), see asterisk in Fig. 6), together with 
its instability eigenvalues. The nonlinear evolution proceeds through periodic 
oscillations between two asymmetric states (see dashed lines for t — 0). 

In the bottom row of Fig. 7 we show the evolution of another unstable solution 
and its associated instability spectrum. In this case, we take the Si solution 
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at (/i, C) = (—0.6,0.4), which is in the gap between the stability regions of 
the Si and S 2 solutions. A loop of pitchfork bifurcations occurs in the gaps 
between consecutive S-stability regions, see details in the next section. The 
evolution of this unstable solution amounts to periodic oscillations between 
the original Si soliton (the solid line at t — 0) and its translation by one site 
(the dashed line at t = 0)). Therefore, this solution is a time-periodic one, 
being close to a heteroclinic connection linking the two unstable solutions, Si 
and its translation. 
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Fig. 8. Homoclinic tangles generating the solitons in Fig. 11. a) the point 
(fJ>,C) = (—0.6,0.1) belongs to the bistability region (see asterisk in Fig. 
6) that generates two solitons (triangles, Ti, and squares, Si), b) the point 
(/i, C) = (—0.7,0.1855) is located at the bifurcation curve (see the circle in Fig. 
6), hence the two solutions collapse into one (see the right panel in Fig. 11). 



4 Bifurcations of the discrete solitons 



In this section we aim to explain how the T\ solutions disappear and how 
the evolution of the S solutions leads to stability changes and creation of 
pairs of stable asymmetric solitons. It is straightforward to understand the 
simultaneous disappearance of the S and T-type solutions as C increases. The 
boundary on which this happens corresponds to the left side of the tent-like 
stability area of the T solutions in Fig. 6. The Sk and Tk solutions collide on 
this boundary and disappear in a saddle-node bifurcation. This bifurcation can 
be easily followed using the homoclinic-tangle approach (see Fig. 8). Before 
the bifurcation occurs, i.e., below the boundary (see the asterisk in Fig. 6), 
two intersections, identified by squares and triangles in Fig. 8(a), generate 
the solitons of the types Si and Ti, respectively, which are are shown in the 
left panels of Fig. 11. In contrast, when the bifurcation curve is approached, 
see Fig. 8(b), the stable and unstable manifolds barely intersect, and the two 
solitons are nearly identical. Exactly at the bifurcation, the manifolds touch 
tangentially, and only one solution is generated (i.e., the Si and Ti solitons 
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are identical at this point). After this saddle-node bifurcation, these solutions 
do not exist anymore. 
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Fig. 9. Different types of bifurcations around the cusp of the T\ region. A sad- 
dle-node bifurcation corresponding to panel (e) annihilates the pair of the Si and 
T\ solutions (the same bifurcation was depicted in Fig. 8). Another saddle-node 
bifurcation corresponds to panel (g) and is responsible for the destruction of the T\ 
solution. Going from (d) to (a), there is a pitchfork bifurcation which is responsible 
for the loss of stability of the symmetric soliton and the creation of a pair of stable 
asymmetric ones (see further details in Fig. 10). 

Another family of noteworthy bifurcations points corresponds to cuspidal 
points of the T-stability regions. Three types of bifurcations occur near these 
points: (A) the saddle-node bifurcation described above, at which the S and T 
solitons collide and disappear; (B) another saddle-node bifurcation, at which 
the T solutions disappear; and (C) a bifurcation at which the S solitons lose 
their stability. These bifurcations are depicted in Fig. 9 as follows: (A) corre- 
sponds to going through panels (f)— >(e)— »(b), (B) is displayed by a chain of 
panels (f)— >(g)— >(d), and (C) corresponds to going through (d)— >(a)— »(b). 

Bifurcation (C) deserves more attention. In Fig. 10 we depict in more detail 
the bifurcation scenario corresponding to a route from the stability region of 
Si to its S2 counterpart, which passes through a gap where no S solution is 
stable. Panels (a) through (g) show the homoclinic tangles as C increases, 
for fixed fi = —0.55. Parameter values for each of these panels are indicated 
in panel (h) and its magnification, panel (i). From this figure, the existence 
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Fig. 10. The pitchfork bifurcation leading to asymmetric soliton solutions. This 
series of plots was obtained by keeping fi = —0.55 and varying C from the Si to 
S2 regions, as shown in insets (h) and (i). The series includes the following panels: 
(a) the Si region before the bifurcation; (b) at the pitchfork bifurcation; (c), (d) 
and (e) in the region where asymmetric solitons exist; (f) is a reverse pitchfork that 
destroys the asymmetric solutions; and (g) the S2 region after the latter bifurcation. 
The square dot represents the same homoclinic orbit throughout the bifurcation, 
while the triangle and circle in (d) correspond to the asymmetric solitons from Fig. 
3. 

of a pitchfork loop — a supercritical pitchfork shortly followed by a reverse 
subcritical one — is evident. In panel (a) we depict by a square point one of 
the homoclinic intersections, which gives rise to solution Si. As C increases, 
the supercritical pitchfork bifurcation is reached, panel (b), that gives rise to 
a pair of extra homoclinic intersections (see panel (c) and triangle and circle 
in panel (d)). As C increases further, the reverse subcritical pitchfork occurs, 
eliminating the extra pair of intersections. 

Quite interesting are the shape and stability of solitons that this extra pair 
generates in the narrow gap between the Sk and Sk+i stability regions (which 
coincides with the region of the pitchfork loop), where these new solutions ex- 
ist. They represent asymmetric solitons of the CQ DNLS equation, which are 
depicted in Fig. 3(b) by means of triangles and circles. As might be expected, 
the stability is swapped by the pitchfork bifurcations. Indeed, the stable S 
solitons are unstable inside the pitchfork loop, where the new asymmetric so- 
lutions are stable. An example of unstable evolution of an S soliton inside 
the pitchfork loop is shown in the bottom row of Fig. 7. We stress that such 
pitchfork loops, and the respective asymmetric solitons, do not exist in the cu- 
bic DNLS equation (asymmetric solitons were also found in a DNLS equation 
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with a mixture of cubic onsite and intersite nonlinearities [19]). 

Further inspection of the pitchfork loop demonstrates that, if the Si site- 
centered soliton looses its stability through the direct bifurcation, then the 
reverse bifurcation, which closes the loop, stabilizes an S2 bond- centered soli- 
ton. At the next pitchfork loop (in the gap between stability of S2 and S3), 
the latter bond-centered solution gets destabilized and, after the loop, an S3 
site-centered soliton gains its stability. As C is increased further, this process 
repeats itself, creating stability bands for the Sk solitons for larger k. In fact, 
all the solutions Sk represent only two distinct soliton families emanating from 
the principal homoclinic intersection: of site-centered solutions, and of bond- 
centered ones, the pitchfork loops conducting stability swaps between the two 
families. For given Sk, the number of the pitchfork loops passed by the soliton, 
while developing from the anti-continuum (C = 0) limit, is A; — 1. 

The existence of stable bond-centered solitons is a noteworthy feature by itself, 
as in the usual DNLS model with the cubic nonlinearity only site-centered 
states are stable (in the DNLS equation with the saturable onsite nonlinearity, 
stable bond-centered solitons were found too [20]). An interesting issue which 
still has to be addressed is whether there exists an accumulation curve for 
the pitchfork loops, or the loops continue to appear as one approaches the 
continuum limit, C = 00. 



5 The variational approximation 

The variational approximation (VA) can be used to describe the shape of 
stationary soliton solutions in an analytical form. The method will not only 
be successful in approximating the shape of the most fundamental solitons, 
but will also predict the saddle-node bifurcation where the Ti and Si solutions 
collide and disappear. 

The only tractable ansatz for the VA in the discrete models is one based on 
an exponential cusp, that was applied in Ref. [25] to solitons in the above- 
mentioned DNLS equation with an arbitrary power nonlinearity : 



with real positive constants A and a. We identify the value of a without the 
resort to the VA proper, but rather from the substitution of ansatz (8) in Eq. 
(5) linearized for the decaying tail of the soliton, which yields 



u n = Ae 



a\n\ 
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(recall a = 2 — (i/O). Thus, a necessary (but not sufficient) condition for the 
existence of any soliton solution is a > 2, or // < (since we set C > 0, 
it will then guarantee that Eq. (9) yields a real positive a). We stress that 
the latter soliton-existence condition is an exact one, as it follows from the 
straightforward consideration of the exponentially localized tails of the soliton, 
and does not exploit any approximation. 

Now we invoke the VA proper, treating the amplitude A in ansatz (8) as 
a variational parameter, while a was already fixed as per Eq. (9) (i.e., by 
the condition that the ansatz must match the correct asymptotic form of the 
soliton solution). 




-4 -2 2 4 -4 -2 2 4 

n n 




-10 10 -10 10 

n n 



Fig. 11. Soliton solutions of Eq. (5) for (fi,C) = (-0.6,0.1) (left) and 
(/i,C) = (—0.7,0.1845) (right) are shown on straight (top) and logarithmic (bot- 
tom) scales. Points and lines correspond, respectively, to numerically exact solutions 
and ones predicted by the variational approximation. The right panels correspond 
to a point (fi, C) close to the collision of the soliton solutions of the S\ and T\ types 
(see the circle in Fig. 6). 

The stationary equation (5) can be derived from the Lagrangian 



E 



B 



Q 



—u n - —u n - C{u n+ i 



u, 



(10) 



The substitution of the ansatz into the Lagrangian and explicit calculation of 
the sum lead to the following effective Lagrangian, which is the main ingredient 
of the VA [30]: 

A 6 

L eS (A) = fiA 2 coth(a) + A 4 coth(2a) coth(3a) - 2 tanh(a/2)CA 2 . (11) 

3 



16 



Then, the Euler-Lagrange equation, dL e s{A) / dA = 0, yields a quadratic equa- 
tion for A 2 , 

-coth(3a) (A 2 ) 2 + 2coth(2a) (A 2 ) +/icoth(a) - 2 tanh(a/2)C = 0. (12) 




0-1 0-1 



Fig. 12. Amplitude A (left) and norm M (right) for the two soliton solutions pre- 
dicted by the VA (the scale for the norm is logarithmic). The top and bottom 
branches correspond to the tall and short solitons, respectively. The dark solid line 
represent the bifurcation curve where the two solutions are predicted to coalesce 
and disappear. 

Equation (12) indicates that there may be two different solutions. The top 
panel of Fig. 12 depicts the amplitude associated with these two solutions 
predicted by the VA. It is clear that the two solutions coalesce and disappear 
at the bifurcation curve (solid line). In Fig. 11 we display examples of these 
two variational solutions, together with the ones obtained numerically from 
Eq. (5) through Newton iterations. It is observed that the match between 
the numerical and variational solutions is extremely good for small values of 
the coupling parameter C. For larger values of C, the solution tends to its 
continuum (smooth) analogue where the cusp-shaped ansatz (8) is obviously 
irrelevant. 

The VA predicts the bistable solutions of the type (8) inside the (//, C)- 
parameter regions where the quadratic equation (12) admits two distinct pos- 
itive roots. The comparison with numerical results in Fig. 11 immediately 
shows that the two different variational solutions exactly correspond to the S\ 
and Ti solitons, as they were defined above. Therefore, the curve where Eq. 
(12) has a double root represents the saddle-node collision of Si and T ± . This 
bifurcation curve predicted by the VA is depicted by the dashed-dotted curve 
labeled VAl in Fig. 6. It is quite remarkable that the approximation based on 
the simple ansatz (8) is able to capture the saddle-node bifurcation so well for 
small C. 

It is possible to refine the VA approach by allowing a more general ansatz of 



17 



the form 

4 e -«M if \n\ > 0, 



(13) 

13 A if n = 0, 



where we introduce a new free parameter (5. This new ansatz is able to predict 
the existence of the solitons of two distinct types, which can me immediately 
identified with short (S) and tall (T) solitons, that were described in detail 
above. 

The new variational ansatz (13), with two independent variational parameters 
(A, (3), is able to very accurately predict the whole existence region of the Ti 
solution, as well as the saddle-node bifurcation which creates the S3 solution. 
The boundary of the £3 solitons predicted by the improved VA is depicted by 
the dash-dotted line VA2 in Fig. 6. As seen from the figure, the new VA gives 
an extremely good approximation for the boundary, up to C < 1. The results 
for the existence region of the Ti soliton predicted by this VA are not depicted 
in Fig. 6 because they exactly coincide (up to the resolution of the figure) with 
both, left and right, numerically found boundaries for the Ti soliton, see the 
darker shaded area in Fig. 6 (the simple VA, based on ansatz (8), which gives 
rise to curve VAl, only captured the left boundary). 

The improved VA with ansatz (13) is also able to pick up solutions with a 
dip at the central site, corresponding to /3 < e~ a in the ansatz, i.e., a bound 
state of two solitons, or multibreather [29]). Thus, the modified VA is able to 
describe the respective bifurcation curves too (results not presented here). 

A well-known approach to predicting the stability for soliton families is based 
on the Vakhitov-Kolokolov (VK) criterion [31]. The VK criterion states that 
a solution family, parameterized by the frequency /i, as in Eq. (4), may be 
stable if dM/d/i < 0, and is definitely unstable otherwise, where M is the 
soliton's norm defined as per Eq. (3). The simplest variational ansatz (8) 
yields M = A 2 coth(a) [note that both a and A depend on (/i, C) through 
Eqs. (9) and (12)]. The norm given by the latter expression is depicted, as a 
function of (//, C), in the right panel of Fig. 12. It is clear from the figure that 
dM/d/j, > for the Ti soliton, and dM/d/j, < for Si (these results can be 
proven analytically, within the framework of the VA). Thus, the VK criterion 
suggests that the Ti soliton is unstable, while its Si counterpart may be stable. 
Comparing this conclusion with the numerical results for the stability reported 
above we conclude that the VK criterion does not apply to our model. In fact, 
similar conclusions for the failure of the VK criterion where obtained in the 
continuous CQ NLS equation with an external potential, that might be both 
a single rectangular potential well [17], and a periodic lattice of rectangular 
wells (the above-mentioned Kronig-Penney potential) [16]. 

A natural question that arises is the existence of multistable solutions in 
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Fig. 13. Two-dimensional coexisting soliton solutions for (fi,C) = (—0.6,0.07). The 
top row depicts three coexisting stable soliton solutions while the bottom row rep- 
resents three coexisting unstable soliton solutions. All these solutions were obtained 
by Newton iteration on the steady state equation with initial seeds predicted by the 
two-dimensional variational approximation with ansatz (14). 



higher dimensional lattices. The higher dimensional equivalent of (1) is ob- 
tained by replacing the one- dimensional index n and the discrete Laplacian by 
their higher- dimensional analogues. Unfortunately, the homoclinic approach is 
not applicable in higher dimensional lattices. Nonetheless, the variational ap- 
proach presented here is still amenable in the higher dimensional case. Namely, 
a natural ansatz equivalent to (13) but now in two dimensions (2D) becomes 

f Ae ~a(\n\ + \m\) jf | n | + | m | > Q) 
U nm = \ (14) 

yf3A if n = m = 0, 

where now we have two indexes (m, n) along the two spatial dimensions, and 
the decay a is the same as for the one-dimensional case [cf. Eq. (9)]. Pre- 
liminary results indicate the coexistence of stable soliton solutions in the 2D 
model. As an example we depict in Fig. 13 six coexisting solutions for the 
same parameter values. The top row in the figure depicts three stable coex- 
isting solutions. It is worth mentioning that the variational approach for the 
2D case is also able to capture other solutions such as, stable and unstable, 
multi-humped profiles (see bottom row of Fig. 13). A detailed stability anal- 
ysis for higher dimensional (2D and 3D) soliton solutions falls outside of the 
scope of the present manuscript and will be presented elsewhere. 
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6 Conclusions and extensions 



The objective of this work was to introduce the discrete nonlinear Schrodinger 
(DNLS) equation with the competing cubic-quintic (CQ) nonlinearities. Be- 
sides being a new dynamical model, this model may also apply to the descrip- 
tion of an array of optical waveguides with the intrinsic CQ nonlinearity. 

We have studied the multistability of discrete solitons in the model, looking for 
homoclinic solutions to the respective stationary discrete equation. Regions of 
the existence and stability of single-humped soliton solutions were identified by 
using a numerical continuation method based on Newton-type iterations. The 
stability of the various types of the solitons was investigated through numer- 
ical evaluation of eigenvalues for small perturbations. The resulting stability 
diagram suggests the existence of an infinite family of branches of stable soli- 
tons of the aforementioned type Tk, for sufficiently small values of the coupling 
constant C . As C increases, these solutions get destroyed through saddle-node 
bifurcations. We have also identified another type of discrete soliton solutions, 
Sk, with Si surviving in the continuum limit, C — > oo. The stability of solu- 
tions of the latter type (S) changes, with the increase of C, through a series 
of small pitchfork loops (each opens with a supercritical pitchfork, which is 
shortly followed by a reverse supercritical one). Inside the loops the symmetric 
solitons Sk loose their stability, while a pair of stable asymmetric solutions is 
created. The latter solutions have no counterpart in the cubic DNLS equation, 
nor in the continuum CQ NLS equation. Finally, using the variational approx- 
imation (VA), we were able to approximate the main branches of the solutions 
and their bifurcations for small C. We have also proposed an improved ver- 
sion of the VA, with two free parameters rather than one, that drastically 
(qualitatively) upgrades the accuracy of the VA, making it possible to predict 
simultaneously the most fundamental solitons and some higher-order ones. 

An interesting extension of this work would be to thoroughly investigate lo- 
calized solutions in the two-dimensional version of the CQ DNLS equation, 
which may be realized, for instance, as a bunch of the corresponding nonlin- 
ear optical waveguides (cf. the recently reported experimental realization for 
linear waveguides [32]). In particular, the existence and stability of asymmet- 
ric solitons in the two-dimensional model would be an issue of great interest. 
Another promising direction for further considerations may be the study of 
kink solutions, corresponding to heteroclinic trajectories generated by the in- 
tersection of the stable and unstable manifolds of different fixed points of the 
map (6). 

As concerns the underlying mathematical theory, an intriguing issue is to 
understand why the Vakhitov-Kolokolov criterion was found to consistently 
fail in the discrete CQ NLS equation and its counterpart with a periodic 
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potential in the continuum case [16,17]. 
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